#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <list>
#include <cmath>

#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 2>::tiny floatTinyArray_t;
typedef TunaArray<double, 2>::huge ScalarField2D;
typedef TunaArray<double, 1>::huge ScalarField1D;

void writeToFile_DX(list<floatTinyArray_t> &, int, string);
void writeToFile_GNU(list<floatTinyArray_t> &, int, string);
void writeToFile_VTK(const ScalarField2D&, std::string, std::string);
void readFromFile_DX_2D(ScalarField2D&, ScalarField2D&, 
			int, string, int w = 5); 

timer crono;

int main () 
{
    double length_x, length_y, dt, left_extreme, right_extreme, tol, dist;
    double s0, cD;
    int num_nodes_x, num_nodes_y, Npoints, frec, Nstart, Nend, jump, Nxf, Nyf;

    double etime = 0;
// -----> Reading data from the input file
    std::ifstream input_cin ("input01");
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> left_extreme
	      >> right_extreme
	      >> Npoints
	      >> Nstart
	      >> Nend
	      >> jump
	      >> dt
	      >> s0
	      >> cD
	      >> tol
	      >> Nxf
	      >> Nyf
	      >> frec;
    input_cin.close();

    /*
    cout << "\n +-----+ \n" 
	 << length_x << "\t" << length_y << "\t" << num_nodes_x << "\t" << num_nodes_y << "\n"
	 << left_extreme << "\t" << right_extreme << "\t" << Npoints << "\n"
	 << Nstart << "\t" << Nend << "\t"<< dt << "\n"
	 << tol << "\t" << frec 
	 << "\n +-----+ \n";
	 /**/

// -----> Mesh for velocity
    StructuredMesh< Uniform<double, 2> > mesh(length_x, num_nodes_x, 
					      length_y, num_nodes_y);
    floatTinyArray_t deltas = mesh.getDeltas();
    mesh.print();
    
// -----> List for points storage
    list<floatTinyArray_t> linec;
    list<floatTinyArray_t>::iterator iter;
    floatTinyArray_t aux; 
    floatTinyArray_t x_i_n;     // Coordinates of x(i) at step n
    floatTinyArray_t x_i_np1;   // Coordinates of x(i) at step n+1
    floatTinyArray_t x_im1_n;   // Coordinates of x(i-1) at step n
    floatTinyArray_t x_im1_np1; // Coordinates of x(i-1) at step n+1
    double deltaL0 = (right_extreme - left_extreme) / (Npoints - 1);
    //    double s0_deltaL0 = s0 * deltaL0;

    double line_length = 0.0;
    ofstream flen("./DataMix2D_01L/lineLength"); 

    cout << "\n +----- Line construction -----+ ... ";

    list<double> s, tao, deltaL;
    list<double>::iterator iS, iT, iD;

    for(int i = 0; i < Npoints; i++) {
	aux(0) = left_extreme + i * deltaL0;
	aux(1) = 0.5;
	linec.push_back(aux);

	// We need an S and a Tao for each point (interval)
	s.push_back(s0);
	tao.push_back(0.0);
	deltaL.push_back(deltaL0);
    }

    cout << "\n Sizes: line = " << linec.size() 
	 << "; s = " << s.size() 
	 << "; tao = " << tao.size()
	 << "; deltaL = " << deltaL.size() << "\n" << flush;

    cout << " +----- DONE -----+\n";

    writeToFile_DX(linec, 0, "./DataMix2D_01L/line.");

    cout << "\n +--------------------------+"
	 << "\n +      Begining track      |"
	 << "\n +--------------------------+\n\n";

// -----> Velocity field at time t
    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());

// -----> Velocity field at time t + dt
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());

    ParticleTracking< RungeKutta4< Linear<double, 2> > > particle;

    ScalarField2D c(Nxf, Nyf);
    double sigmax, sigmay, normSig, coordx, coordy, midx, midy;
    double aux1, aux2, auxS, auxT, auxD;
    double dxf = length_x / (Nxf - 1);
    double dyf = length_y / (Nyf - 1);

    cout << " +----- Time series -----+\n"
	 << " | Nstar = " << Nstart << "\t Nend = " << Nend 
	 << "\t jump = " << jump 
	 << "\n +-----------------------+ ";
// -----> Main loop

    int pausa, contador = 0;
    for(int step = Nstart; step < Nend; step += jump) {
      

      cout << "\n\n ----> Step = " << step << "\t Npoints = " << linec.size() << flush;
      cout << "\n Sizes: line = " << linec.size() 
	   << "; s = " << s.size() 
	   << "; tao = " << tao.size()
	   << "; deltaL = " << deltaL.size() << "\n" << flush;
      
      cout << "\n   +----- Tracking ... " << flush;
      crono.tic();
// -----> Read the velocity at time t from a file
      readFromFile_DX_2D(u, v, step, "/home/luiggi/DataMix2D_01/velc.");
//      InOut::readFromFile_1(u, v, step, "../Velocity/Data01/velc.");
// -----> Read the velocity at time t + dt from a file
      readFromFile_DX_2D(un, vn, step + jump, "/home/luiggi/DataMix2D_01/velc.");
//      InOut::readFromFile_1(un, vn, step + jump, "../Velocity/Data01/velc.");

// -----> Track every particle on the line

      iter = linec.begin(); // Pointer to initial point in the line x(i-1)^n
      x_im1_n = *iter;      // Contents of initial point in the line x(i-1)^n

// Track first point to the next position x(i-1)^(n+1) 
      *iter = particle.track(*iter, deltas, u, v, un, vn, dt);
      x_im1_np1 = *iter;
      iter++;  // Advance iter to the second point in the line

      line_length = 0.0;

      contador = 0;
      
      iD = deltaL.begin(); iS = s.begin(); iT = tao.begin();
      iD++; iS++; iT++;

      for( ; //iD = deltaL.begin(), iS = s.begin(), iT = tao.begin();
	   iter != linec.end(); 
	   ++iter, ++iS, ++iT, ++iD)  {

	x_i_n = *iter; // Current point  x(i)^n
	x_i_np1 = particle.track(*iter, deltas, u, v, un, vn, dt); 

/// Calculate the distance between consecutive tracked points 
/// x(i-1)^(n+1) and x(i)^(n+1)
	dist = sqrt( ( x_i_np1(0) - x_im1_np1(0) ) * 
		     ( x_i_np1(0) - x_im1_np1(0) ) +
		     ( x_i_np1(1) - x_im1_np1(1) ) * 
		     ( x_i_np1(1) - x_im1_np1(1) ) );

/// if dist exceeds the tolerance, adds a point between points: 
/// x(i)^n and x(i-1)^n

	if ( !(step % 10) ) {

	  
	  aux(0) = ( x_i_n(0) + x_im1_n(0) ) * 0.5;
	  aux(1) = ( x_i_n(1) + x_im1_n(1) ) * 0.5;
	  linec.insert(iter,aux); // Insert one point, before the actual position of iter
	  iter--; iter--; // Rewind the iterator 2 positions, (the for loop will advance one pos.)

	  deltaL.insert(iD, dist * 0.5); // Insert one delta, before the actual position of iD
	  *iD = dist *0.5; // Update the delta in the actual position of iD
	  iD--;iD--; // Rewind iD one position, to the previously inserted delta

	  s.insert(iS, *iS); // Striation thickness will keep its value for the new two intervals.
	  iS--;iS--; // Rewind iD one position, to the previously inserted striation thickness

	  tao.insert(iT, *iT); // Dimensionless time will keep its value for the new two intervals.
	  iT--;iT--; // Rewind iD one position, to the previously inserted tao.
	  /* *
	  if (step == 1020 || step == 990) {
	    cout << "\n count = " << contador++ << " JAJAJA \n" << flush;
	    cout << "\n Sizes: line = " << linec.size() 
		 << "; s = " << s.size() 
		 << "; tao = " << tao.size()
		 << "; deltaL = " << deltaL.size() << "\n" << flush;
	    cout << " s = " << *iS << "; tao = " << *iT << "; deltaL = " << *iD << "\n" << flush;
	  }
	  /* */
	} else { // If dist does not exceeds the tolerance, 
	  /* *
	  if (step == 1020 || step == 990) {
	    cout << "\n count = " << contador++ << " XXXXX \n" << flush;
	    cout << "\n Sizes: line = " << linec.size() 
		 << "; s = " << s.size() 
		 << "; tao = " << tao.size()
		 << "; deltaL = " << deltaL.size() << "\n" << flush;
	    cout << " s = " << *iS << "; tao = " << *iT << "; deltaL = " << *iD << "\n" << flush;
	  }
	  /* */
	  x_im1_n   = x_i_n;   // Re-state the value of x(i-1)^n
	  x_im1_np1 = x_i_np1; // Re-state the value of x(i-1)^(n+1)
	  *iter     = x_i_np1; // Modify the new position of x(i) in the list ( x(i)^(n+1) )
	  
	  *iS =  (*iS) * (*iD) / dist; // Striation thickness: s_i = s0 * delta_0 / delta_i
	  *iT = *iT + dt * cD / ( (*iS) * (*iS) ); // Calc. dimensionless time
	  *iD = dist; // Update the delta_i

	  line_length += dist; // Accumulate the length of the strip
	}
      }
      cout << " Elapsed time = " << crono.toc() << flush ;

      flen << step << "\t" << line_length << endl;
 
      if( !(step % frec) ) {
	cout << "\n   +----- Writing file ... " << flush ;
	crono.tic();
	writeToFile_DX(linec, step, "./DataMix2D_01L/line.");
	cout << " Elapsed time = " << crono.toc() << flush ;
      }
      cout << "\n   +----- Diffusion calculation ... " << flush;

      /* *
// ------> 1) Calc. striation thickness and dimensionless time
      crono.tic();

      iter = linec.begin();  
      iS = s.begin();
      iT = tao.begin();
      x_im1_n = *iter; iter++; 
      
      deltaL0 = 0.0;

      for( ; iter != linec.end(); ++iter, ++iS, ++iT) {
	x_i_n = *iter;
	*iS = s0_deltaL0 / sqrt( (x_i_n(0) - x_im1_n(0)) * (x_i_n(0) - x_im1_n(0)) +
				 (x_i_n(1) - x_im1_n(1)) * (x_i_n(1) - x_im1_n(1)) );
	*iT = *iT + dt * cD / ( (*iS) * (*iS) );
	deltaL0 += *iS * sqrt(1 + 4 * (*iT));
	x_im1_n = *iter;
      }    
      
      deltaL0 /= linec.size();
      /* */
// ------> Reconstruct the concentration field on a fine grid 
      c = 0;
      for(int i = 0; i < Nxf; ++i)
	for(int j = 0; j < Nyf; ++j) {
	  iter = linec.begin();  
	  x_im1_n = *iter; 
	  iter++;
	  for( iS = s.begin(), iT = tao.begin(), iD = deltaL.begin(); 
	       iter != linec.end(); 
	       ++iter, ++iS, ++iT, ++iD) {
	  
	    x_i_n = *iter;
	    sigmax = x_i_n(0) - x_im1_n(0); 
	    sigmay = x_i_n(1) - x_im1_n(1);
	    normSig = sqrt ( sigmax * sigmax + sigmay * sigmay);
	    sigmax /= normSig;
	    sigmay /= normSig;
	    coordx = dxf * i;
	    coordy = dyf * j;
	    midx = 0.5 * ( x_i_n(0) + x_im1_n(0) ); 
	    midy = 0.5 * ( x_i_n(1) + x_im1_n(1) );	  
	    aux1 = (coordx - midx) * sigmax + (coordy - midy) * sigmay;
	    aux2 = (coordx - midx) * (-sigmay) + (coordy - midy) * sigmax;
	    aux1 = aux1 * aux1 / ( (*iD) * (*iD) );
//( (*iS) * (*iS) * (1.0 + 4.0 * (*iT)) );
//(deltaL0 * deltaL0);
	    aux2 = aux2 * aux2 / ( (*iS) * (*iS) * (1.0 + 4.0 * (*iT)) );
	    c(i,j) += exp(-aux1-aux2) / (1.7726 * sqrt( 1.0 + 4.0 * (*iT)));
	    x_im1_n = *iter;
	  }
	}
      cout << " Elapsed time = " << crono.toc() << flush ;
      
      cout << "\n   +----- Writing file ...." << flush;
      crono.tic();
      InOut::writeToFile_DX(c, step, "./DataMix2D_01L/conc.", dxf, dyf); 	    
      //writeToFile_VTK(c, "./DataMix01/conc.vtk", "Concen"); 
      cout << " Elapsed time = " << crono.toc() << flush;

      /* */

    }

    flen.close();
    cout << "\n End ... \n\n";

    return 0;
}

void writeToFile_DX(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    //    cout << "writing " << filename << endl;
    
    file << " object 1 class array type float rank 1 shape 2 items "
	 << Npoints << " data follows\n";
    for(iter = linec.begin(); iter != linec.end(); ++iter)
	file << (*iter)(0) << "   " << (*iter)(1) << "   ";
    file << endl;

    file << "object 2 class array type int rank 1 shape 2 items "
	 << Npoints - 1 << " data follows\n";
    for(int i = 0; i < Npoints - 1; ++i)
	file << i << "   " << i + 1 << "   ";
    file << endl;

    file << "attribute \"element type\" string \"lines\"\n"
	 << "attribute \"ref\" string \"positions\"\n";
    file << "object 3 class array type float rank 0 items "
	 << Npoints << " data follows\n";
    for(int i = 0; i < Npoints; ++i)
	file << 1 << "\t";
    file << endl;
 
    file << "attribute \"dep\" string \"positions\"\n";
    file << "object \"irregular positions regular connections\" class field\n"
	 << "component \"positions\" value 1\n"
	 << "component \"connections\" value 2\n" 
	 << "component \"data\" value 3\n" 
	 << "end\n" ;
    file.close();

}

void writeToFile_GNU(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    cout << "writing " << filename << endl;
    
    for(iter = linec.begin(); iter != linec.end(); ++iter)
	file << (*iter)(0) << "\t" << (*iter)(1) << "\n";
    file << endl;

    file.close();

}

void writeToFile_VTK(const ScalarField2D& data, 
		     std::string name, std::string variable)
{
    int Nvtk_x = data.rows();
    int Nvtk_y = data.cols();
    int Nvtk_z = 1;
    double dx = 1. / (Nvtk_x - 1);
    double dy = 1. / (Nvtk_y - 1);
    double dz = 1.;
    int Ntotal = Nvtk_x * Nvtk_y * Nvtk_z;

    std::ofstream file (name.c_str());        
    file << "# vtk DataFile Version 2.0 \n"
	 << "Lid_driven streamfuntion-vorticity \n"
	 << "ASCII \n"
	 << "DATASET RECTILINEAR_GRID \n" 
	 << "DIMENSIONS " 
	 << Nvtk_x << " " << Nvtk_y << " " << Nvtk_z;

    file << "\nX_COORDINATES " << Nvtk_x << " float \n";
    for(int i = 0; i < Nvtk_x; ++i) {
	file << i * dx << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nY_COORDINATES " << Nvtk_y << " float \n"; 
    for(int i = 0; i < Nvtk_y; ++i) {
	file << i * dy << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nZ_COORDINATES " << Nvtk_z << " float \n"; 
    for(int i = 0; i < Nvtk_z; ++i) {
	file << i * dz << " ";
	if ( !(i % 6) ) file << "\n";
    }
	
    file << "\nPOINT_DATA " << Ntotal
	 << "\nSCALARS " << variable << " float"
	 << "\nLOOKUP_TABLE default \n";
    for(int k = 1; k <= Nvtk_z ; ++k) {
	for(int j = 1; j <= Nvtk_y ; ++j) {
	    for(int i = 1; i <= Nvtk_x; ++i) {
		file << data(i,j) << " ";
		if ( !(i % 6) ) file << "\n";
	    }
	}
    }
        
    file.close();

}

void readFromFile_DX_2D(ScalarField2D& scf1, ScalarField2D& scf2, 
			int number, string filename, int w) 
{   
  //  int w = 5;
  ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;
    ifstream file (filename.c_str());

    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    } else {
      //      cout << "\n Opening file \" " << filename << "\"" << endl;
    }
    
    char line[256];
/// 
/// by the time, just jump the first 8 lines
///	
    for(int ii = 0; ii < 8; ++ii) {
	file.getline(line,256);
	//	cout << " ii = " << ii << "\t" << line << endl;
    }

    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);
///
/// Here, scf1 and scf2 are supposed to be of type ScalarField2D
///
    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j)
	  file >> scf1(i,j) >> scf2(i,j) ;
    file.close();
}


